home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / a_utils / expanded.lha / expanded / src / Map.C < prev    next >
C/C++ Source or Header  |  1992-03-19  |  51KB  |  1,794 lines
  1. //
  2. // Linear-Affine-Projective Geometry Package
  3. //
  4. // Map.C
  5. //
  6. // $Header$
  7. //
  8. // William J.R. Longabaugh 
  9. // University of Washington
  10. //
  11. // Implementation of the linear-affine-projective geometry
  12. // package described in William J.R. Longabaugh, "An Expanded
  13. // System for Coordinate-Free Geometric Programming", Master's 
  14. // thesis, University of Washington, 1992.
  15. //
  16. // Copyright (c) 1992, William J.R. Longabaugh
  17. //   Copying, use and development for non-commercial purposes permitted.
  18. //   All rights for commercial use reserved.
  19. //   This software is unsupported and without warranty; it is
  20. //   provided "as is".
  21. //
  22. // ***********************************************************************
  23.  
  24. #include <math.h>
  25. #include "Lap.h"
  26.  
  27. // ***********************************************************************
  28. // ***********************************************************************
  29. //
  30. // Map Class
  31. //
  32. // ***********************************************************************
  33. // ***********************************************************************
  34. //
  35. // Private/protected member functions
  36. //
  37. // ***********************************************************************
  38. //
  39. // This routine is used by map constructors to convert lists
  40. // of arbitrary geometric objects into objects of the desired type
  41. // in the desired space.  A matrix containing std. coords of the
  42. // objects is built:
  43. //
  44.  
  45. Boolean Map::ConvertList(GeObList& v, GeObType targ,
  46.              SubSet& dest, Matrix* temp)
  47. {
  48.   Space destspace = dest.EmbeddingSpace();
  49.  
  50. // Determine if the matrix reps. need to be normalized:
  51.  
  52.   Boolean norm = 
  53.     ((destspace.Holds() == PROJ_SPACE) && (dest.Holds() == AFFINE_SUBSET));
  54.  
  55.   for (int i = 0; i < v.Length(); i++) {
  56.     GeOb tempgeob = v[i].MapTo(targ);
  57.     if (tempgeob.SpaceOf() != destspace) {
  58.       return (FALSE);
  59.     } 
  60.     if (!dest.IsIn(tempgeob)) {
  61.       return (FALSE);
  62.     }
  63.     if (norm) {
  64.       tempgeob = dest.Normalize(tempgeob);
  65.     }
  66.     (*temp)[i] = tempgeob.MatrixRep()[0];
  67.   }
  68.   return (TRUE);
  69. }
  70.  
  71. // ***********************************************************************
  72. //
  73. // Apply associated linear map (without explicitly building it):
  74. //
  75.  
  76. GeOb Map::ApplyAL(GeOb& v)
  77. {
  78.   static char* sig = "GeOb Map::ApplyAL(void)";
  79.  
  80.   GeOb temp = v.MapTo(AFF_VECTOR);
  81.   if (temp.SpaceOf() != domain.TangentSub().EmbeddingSpace()) {
  82.     errh.ErrorExit(sig, "Object cannot be mapped to domain space", *this, v);
  83.   } 
  84.   if (!domain.TangentSub().IsIn(temp)) {
  85.     errh.ErrorExit(sig, "Mapped object not in specified domain subset",
  86.                    *this, temp);
  87.   }
  88.  
  89. // The range depends on the character of the affine map.  If it is
  90. // to an affine subset of an affine space, the range is the
  91. // corresponding subspace of the tangent space.  If it is to an
  92. // affine subset of a vector space, the range is a parallel vector
  93. // subspace of the vector space.  If it is to a vector subset,
  94. // this range is the same vector subset.
  95.  
  96.   if (range.EmbeddingSpace().Holds() == AFF_SPACE) {
  97.     Matrix alt = temp.MatrixRep() *
  98.              domain.EmbeddingSpace().HoistTangent() * t;
  99.     Matrix holdmatr = range.EmbeddingSpace().HoistTangent();
  100.     GeOb retval(AFF_VECTOR, range.TangentSub().EmbeddingSpace(),
  101.             alt * Transpose(holdmatr));
  102.     return retval;  
  103.   } else {
  104.     Matrix alt = temp.MatrixRep() *
  105.              domain.EmbeddingSpace().HoistTangent() * t;
  106.     GeOb retval(VECTOR, range.EmbeddingSpace(), alt);
  107.     return retval;  
  108.   }
  109. }
  110.  
  111. // ***********************************************************************
  112. //
  113. // Given all the information, builds a new Map
  114. //
  115.  
  116. Map::Map(MapType ty, Boolean i, SubSet& r,
  117.      SubSet& d, Matrix& tm, Matrix& it, Boolean isdf,
  118.      Boolean hal, GeObType dt, GeObType rt) :  range(r), domain(d),
  119.                                t(tm), invt(it)
  120. {
  121.   invertible = i;
  122.   type = ANY_MAP;
  123.   holds = ty;
  124.   hasAL = hal;
  125.   isDefault = isdf;
  126.   domtype = dt;
  127.   rettype = rt;
  128. }
  129.  
  130. // ***********************************************************************
  131. //
  132. // Dumps out all the data for a map:
  133. //
  134.  
  135. void Map::data_out(ostream& c, int indent, char* name) 
  136. {
  137.   char *ibloc = new char[indent + 1];
  138.   for (int i = 0; i < indent; i++) {
  139.     *(ibloc + i) = ' ';
  140.   }
  141.   *(ibloc + indent) = '\0';
  142.  
  143.   c << ibloc << ast;
  144.   c << ibloc << name;
  145.   c << ibloc << "Type of map: ";
  146.   MapTypeOut(c, type);
  147.   c << "\n";
  148.   c << ibloc << "Currently holds: ";
  149.   MapTypeOut(c, holds);
  150.   c << "\n";
  151.   if (holds != NULL_MAP) { 
  152.     c << ibloc << "Domain subset:\n";
  153.     domain.debug_out(c, indent + ERR_IND);
  154.     c << ibloc << "Range subset:\n";
  155.     range.debug_out(c, indent + ERR_IND);
  156.     c << ibloc << "Matrix representation of map:\n";
  157.     t.debug_out(c, indent + ERR_IND);
  158.     c << ibloc << "Type of domain argument: ";
  159.     GeObTypeOut(c, domtype);
  160.     c << "\n" << ibloc << "Type of result: ";
  161.     GeObTypeOut(c, rettype);
  162.     if (invertible) { 
  163.       c << "\n" << ibloc << "Matrix representation of inverse:\n";
  164.       invt.debug_out(c, indent + ERR_IND);
  165.     } else {
  166.       c << "\n" << ibloc << "Map is not invertible\n";
  167.     }  
  168.     if (hasAL) { 
  169.       c << ibloc << "Map has an associated linear map\n";
  170.     } else {
  171.       c << ibloc << "Map does not have an associated linear map\n";
  172.     }
  173.     if (isDefault) { 
  174.       c << ibloc << "Map is a default map\n";
  175.     } else {
  176.       c << ibloc << "Map is not a default map\n";
  177.     }
  178.   }  
  179.   c << ibloc << ast;
  180.  
  181.   delete ibloc;
  182.   return;
  183. }
  184.  
  185. // ***********************************************************************
  186. //
  187. // Public member functions
  188. //
  189. // ***********************************************************************
  190. //
  191. // Need to do memberwise initialization, since matrix class members need
  192. // to do memberwise initialization:
  193. //
  194.  
  195. Map::Map(Map &m) : range(m.range), domain(m.domain), t(m.t), invt(m.invt)
  196. {
  197.   invertible = m.invertible;
  198.   type = ANY_MAP;
  199.   holds = m.holds;
  200.   hasAL = m.hasAL;
  201.   isDefault = m.isDefault;
  202.   domtype = m.domtype;
  203.   rettype = m.rettype;
  204. }
  205.  
  206. // ***********************************************************************
  207. //
  208. // Need to do memberwise assignment, since Matrix class members need to
  209. // do memberwise assignment:
  210. //
  211.  
  212. Map& Map::operator=(Map &m)
  213. {
  214. //
  215. // This weird checking is required to bypass the apparent inheritance of
  216. // assignment under g++ 1.37:
  217. //
  218.   if ((type != ANY_MAP) && (type != m.holds) && (m.holds != NULL_MAP)) {
  219.     errh.ErrorExit("Map& Map::operator=(Map&)",
  220.                    "Illegal assignment attempted", *this, m);
  221.   }
  222.   range = m.range;
  223.   domain = m.domain;
  224.   t = m.t;
  225.   invt = m.invt;
  226.   invertible = m.invertible;
  227.   holds = m.holds;
  228.   hasAL = m.hasAL;
  229.   isDefault = m.isDefault;
  230.   domtype = m.domtype;
  231.   rettype = m.rettype;
  232.   return (*this);
  233. }
  234.  
  235. // ***********************************************************************
  236. //
  237. // Output for debugging:
  238. //
  239.  
  240. void Map::debug_out(ostream &c, int indent) 
  241. {
  242.   static char* name = "Map Object\n";
  243.   this->data_out(c, indent, name);
  244.   return;
  245. }
  246.  
  247. // ***********************************************************************
  248. //
  249. // If the MultiMap has an atomic VSpace or ASpace for a domain, it
  250. // can be converted to a map.  Alternately, if it is a fully evaluated
  251. // map that evaluates to a vector, it can be converted to a map from
  252. // the dual space to the Reals.
  253. //
  254.  
  255. Map::Map(MultiMap &m)
  256. {
  257.   static char* sig = "Map::Map(MultiMap&)";
  258.  
  259.   type = ANY_MAP;
  260.  
  261.   if (m.FullyEvaluated()) {
  262.     *this = Map(GeOb(m));
  263.   } else {
  264.  
  265. //  Check that the domain is an atomic VSpace or ASpace:
  266.  
  267.     if (m.DomainSpace().CPSpaceSize() != 1) {
  268.       errh.ErrorExit(sig, "Domain of multimap must be an atomic space", m);
  269.     }
  270.  
  271.     holds = (m.Holds() == MULTI_LINEAR) ? LIN_MAP : AFF_MAP;
  272.  
  273.     range = m.RangeSpace().FullSet();
  274.     domain = m.DomainSpace().FullSet();
  275.     domtype = m.DomainType();
  276.     rettype = m.RetType();
  277.     isDefault = FALSE;
  278.     hasAL = (domtype == AFF_POINT);
  279.  
  280.     int ddim = m.DomainSpace().MatrixSize();
  281.     int rdim = m.RangeSpace().MatrixSize();
  282.  
  283.     t = Matrix(ddim, rdim);
  284.     for (int i = 0; i < ddim; i++) {
  285.       t[i] = (m.GetStdImage(i))[0];
  286.     }
  287.  
  288.     if ((ddim == rdim) && (fabs(Det(t)) > EPSILON)) {
  289.       invertible = TRUE;
  290.       invt = Inverse(t);
  291.     } else {
  292.       invertible = FALSE;
  293.     }
  294.   }
  295. }
  296.  
  297. // ***********************************************************************
  298. //
  299. // If the GeOb is a vector, it can be converted to a map from the
  300. // dual space to the Reals.
  301. //
  302.  
  303. Map::Map(GeOb& g)
  304. {
  305.   static char* sig = "Map::Map(GeOb&)";
  306.   GeOb temp;
  307.  
  308.   type = ANY_MAP;
  309.  
  310.   if ((g.Holds() == VECTOR) || (g.Holds() == AFF_VECTOR)) {
  311.     temp = g;
  312.   } else if (g.CanMapTo(VECTOR)) {
  313.     temp = g.MapTo(VECTOR);
  314.   } else {
  315.     errh.ErrorExit(sig, "GeOb cannot be mapped to a vector", g);
  316.   }
  317.  
  318.   range = Reals.FullSet();
  319.   domain = temp.SpaceOf().Dual().FullSet();
  320.   t = Transpose(temp.MatrixRep());
  321.  
  322.   if (domain.EmbeddingSpace().Dim() == 1) {
  323.     invt = Inverse(t);
  324.     invertible = TRUE;
  325.   } else {
  326.     invertible = FALSE;
  327.   }
  328.  
  329.   holds = LIN_MAP;
  330.   hasAL = FALSE;
  331.   isDefault = FALSE;
  332.   domtype = domain.EmbeddingSpace().NativeType();
  333.   rettype = VECTOR;
  334. } 
  335.  
  336. // ***********************************************************************
  337. //
  338. // Map application.
  339. //
  340. // It probably makes sense to also provide a special map application operator
  341. // (to the implementing classes ONLY) that assumes its argument is
  342. // the correct type and fills in a result at a location specified
  343. // by an argument pointer.  This would help standard mappings be more
  344. // efficient.
  345. //
  346.  
  347. GeOb Map::operator()(GeOb& v)
  348. {
  349.   static char* sig = "GeOb Map::operator()(GeOb&)";
  350.  
  351.   if (v.CanMapTo(domtype)) {
  352.     GeOb temp = v.MapTo(domtype);
  353.     if (temp.SpaceOf() != domain.EmbeddingSpace()) {
  354.       errh.ErrorExit(sig, "Object cannot be mapped to domain space", *this, v);
  355.     } 
  356.     if (!domain.IsIn(temp)) {
  357.       errh.ErrorExit(sig, "Mapped object not in specified domain subset",
  358.                   *this, temp);
  359.     }
  360.     if (domtype == PROJ_POINT && holds == AFF_MAP) {
  361.       temp = domain.Normalize(temp);
  362.     }
  363.     Matrix retmat = temp.MatrixRep();
  364.     if (!isDefault) {
  365.       retmat = retmat * t;
  366.     }
  367.     GeOb retval(rettype, range.EmbeddingSpace(), retmat);
  368.     return (retval);
  369.   } else if (hasAL && v.CanMapTo(AFF_VECTOR)) {
  370.     return (this->ApplyAL(v));
  371.   } else {
  372.     errh.ErrorExit(sig, "Object cannot be mapped into domain space",
  373.                *this, v);
  374.   } 
  375. }
  376.  
  377. // ***********************************************************************
  378. //
  379. // Since inverses are computed when the map is created, this is cheap. 
  380. //
  381.  
  382. Map Map::Inv(void)
  383. {
  384.   static char* sig = "Map Map::Inv(void)";
  385.  
  386.   if (!invertible) {
  387.     errh.ErrorExit(sig, "Map is not invertible", *this);
  388.   }
  389.   Boolean hal = ((holds == AFF_MAP) &&
  390.          (domtype != PROJ_POINT) &&
  391.          (rettype == AFF_POINT));
  392.   Map retval(holds, invertible, domain, range, invt, t,
  393.          isDefault, hal, rettype, domtype);
  394.   return (retval);
  395. }
  396.  
  397. // ***********************************************************************
  398. //
  399. // Composition of two maps of the same type:
  400. //
  401.  
  402. Map Map::Compose(Map &m)
  403. {
  404.   static char* sig = "Map Map::Compose(Map&)";
  405.   Scalar factor = 1.0;
  406.  
  407.   if (holds == m.Holds()) {
  408.     if (!domain.IsSubset(m.Range())) {
  409.       errh.ErrorExit(sig,
  410.              "Range of first map not a subset of domain of second map",
  411.              *this, m);
  412.     }
  413.  
  414. // Affine maps into/outof projective spaces need a factor to account
  415. // for the need to normalize the offsets of the subsets:
  416.  
  417.     if ((domain.EmbeddingSpace().Holds() == PROJ_SPACE) &&
  418.        (domain.Holds() == AFFINE_SUBSET)) {
  419.       factor = domain.NormVal(m.Range());
  420.     } 
  421.     
  422.     Boolean cinv = invertible && m.Invertible() && m.Range().IsSubset(domain);
  423.     Matrix ct = factor * m.t * t;
  424.     Matrix cinvt;
  425.     if (cinv) {
  426.       cinvt = (1.0 / factor) * invt * m.invt;
  427.     }
  428.     GeObType cdt = m.DomainType();
  429.     GeObType crt = rettype;
  430.     Boolean chal = ((holds == AFF_MAP) &&
  431.             (cdt == AFF_POINT) && (crt != PROJ_POINT)); 
  432.     Boolean cisdf = isDefault && m.isDefault;
  433.     Map retval(holds, cinv, range, m.Domain(), ct, cinvt, 
  434.            cisdf, chal, cdt, crt);
  435.     return (retval);
  436.   } else {
  437.     errh.ErrorExit(sig, "Cannot compose different map types", *this, m);
  438.   }
  439. }
  440.  
  441. // ***********************************************************************
  442. //
  443. // A common operation in graphics programming is to compose affine
  444. // and projective maps to get a projective map.  This operation
  445. // does this composition by automatically casting any affine maps
  446. // to projective maps:
  447. //
  448.  
  449. ProjectiveMap Map::ComposeProj(Map &m)
  450. {
  451.   ProjectiveMap first;
  452.   ProjectiveMap second;
  453.   
  454.   if (holds == PROJ_MAP) {
  455.     first = *this;
  456.   } else {
  457.     first = this->InducedProjective();
  458.   }
  459.  
  460.   if (m.Holds() == PROJ_MAP) {
  461.     second = m;
  462.   } else {
  463.     second = m.InducedProjective();
  464.   }
  465.  
  466.   return (first.Compose(second));
  467. }
  468.  
  469. // ***********************************************************************
  470. //
  471. // Given an affine map between two full affine spaces, we frequently desire
  472. // the corresponding projective map between the neighboring projective
  473. // completion spaces: 
  474. //
  475.  
  476. ProjectiveMap Map::InducedProjective(void)
  477. {
  478.   static char* sig = "ProjectiveMap Map::InducedProjective(void)";
  479.  
  480. // For the current implementation, restrict this operation to affine
  481. // maps between full affine subsets of affine spaces.
  482.  
  483.   if (holds != AFF_MAP) {
  484.     errh.ErrorExit(sig,
  485.     "Can only take induced projective map for an affine map", *this);
  486.   }
  487.  
  488.   if ((range.EmbeddingSpace().Holds() != AFF_SPACE ) ||
  489.       (domain.EmbeddingSpace().Holds() != AFF_SPACE ) ||
  490.       (!range.IsFullSpace()) || (!domain.IsFullSpace())) {
  491.     errh.ErrorExit(sig,
  492.            "Only implemented for maps between full affine spaces",
  493.            *this);
  494.   }
  495.  
  496. // Require that the map be invertible; otherwise, we would need to 
  497. // build a special domain subset for the projective map:
  498.  
  499.   if (!invertible) {
  500.     errh.ErrorExit(sig, "Map must be invertible", *this);
  501.   }
  502.  
  503.   SubSet td = domain.EmbeddingSpace().GetSpace(PROJECT_COMP).FullSet();
  504.   SubSet tr = range.EmbeddingSpace().GetSpace(PROJECT_COMP).FullSet();
  505.   GeObType trt = PROJ_POINT;
  506.   GeObType tdt = PROJ_POINT;
  507.  
  508.   Matrix newt = td.EmbeddingSpace().AffineMapTo(AFFINE).t * t *
  509.                 range.EmbeddingSpace().AffineMapTo(PROJECT_COMP).t;
  510.  
  511.   Matrix newinvt = tr.EmbeddingSpace().AffineMapTo(AFFINE).t * invt *
  512.                    domain.EmbeddingSpace().AffineMapTo(PROJECT_COMP).t;
  513.  
  514.   Map retval(PROJ_MAP, invertible, tr, td, newt, newinvt,
  515.          FALSE, FALSE, tdt, trt);
  516.   return retval;  
  517. }
  518.  
  519. // ***********************************************************************
  520. //
  521. // Given an affine map between two full affine spaces, this gives the
  522. // corresponding linear map between the neighboring linearization spaces: 
  523. //
  524.  
  525. LinearMap Map::InducedLinear(void)
  526. {
  527.   static char* sig = "LinearMap Map::InducedLinear(void)";
  528.  
  529. // For the current implementation, restrict this operation to affine
  530. // maps between full affine subsets of affine spaces.
  531.  
  532.   if (holds != AFF_MAP) {
  533.     errh.ErrorExit(sig,
  534.     "Can only take induced linear map for an affine map", *this);
  535.   }
  536.  
  537.   if ((range.EmbeddingSpace().Holds() != AFF_SPACE ) ||
  538.       (domain.EmbeddingSpace().Holds() != AFF_SPACE ) ||
  539.       (!range.IsFullSpace()) || (!domain.IsFullSpace())) {
  540.     errh.ErrorExit(sig,
  541.            "Only implemented for maps between full affine spaces",
  542.            *this);
  543.   }
  544.  
  545.   SubSet td = domain.EmbeddingSpace().GetSpace(LINEARIZATION).FullSet();
  546.   SubSet tr = range.EmbeddingSpace().GetSpace(LINEARIZATION).FullSet();
  547.   GeObType trt = VECTOR;
  548.   GeObType tdt = VECTOR;
  549.  
  550.   Matrix newt = td.EmbeddingSpace().AffineMapTo(AFFINE).t * t *
  551.                 range.EmbeddingSpace().AffineMapTo(LINEARIZATION).t;
  552.  
  553.   Matrix newinvt;
  554.   if (invertible) {
  555.     newinvt = tr.EmbeddingSpace().AffineMapTo(AFFINE).t * invt *
  556.               domain.EmbeddingSpace().AffineMapTo(LINEARIZATION).t;
  557.   }
  558.  
  559.   Map retval(LIN_MAP, invertible, tr, td, newt, newinvt, 
  560.          FALSE, FALSE, tdt, trt);
  561.   return retval;  
  562. }
  563.  
  564. // ***********************************************************************
  565. //
  566. // The transpose of a map V -> U is a map U' -> V'.  Currently only
  567. // implemented for maps between full spaces.
  568. //
  569.  
  570. LinearMap Map::Trans(void)
  571. {
  572.   static char* sig = "LinearMap Map::Trans(void)";
  573.  
  574. // Make sure this map is linear:
  575.  
  576.   if (holds != LIN_MAP) {
  577.     errh.ErrorExit(sig, "Can only take transpose of a linear map", *this);
  578.   }
  579.  
  580.   if (!range.IsFullSpace() || !domain.IsFullSpace()) {
  581.     errh.ErrorExit(sig,
  582.            "Only implemented for maps between full spaces", *this);
  583.   }
  584.  
  585. // Need to get the dual subspaces for the range and domain,
  586. // and swap them:
  587.  
  588.   SubSet td = range.EmbeddingSpace().Dual().FullSet();
  589.   SubSet tr = domain.EmbeddingSpace().Dual().FullSet();
  590.   GeObType trt = tr.EmbeddingSpace().NativeType();
  591.   GeObType tdt = td.EmbeddingSpace().NativeType();
  592.  
  593.   Map retval(LIN_MAP, invertible, tr, td, Transpose(t),
  594.          Transpose(invt), FALSE, FALSE, tdt, trt);
  595.   return retval;  
  596. }
  597.  
  598. // ***********************************************************************
  599. //
  600. // If this is map A -> X, with A a subset of an affine space and X a
  601. // subset of a linear or affine space, this is the map induced on the
  602. // tangent space A.v: 
  603. //
  604.  
  605. LinearMap Map::AssocLinear(void)
  606. {
  607.   static char* sig = "LinearMap Map::AssocLinear(void)";
  608.  
  609. // First need to make sure that the map has an associated linear map:
  610.  
  611.   if (!hasAL) {
  612.     errh.ErrorExit(sig, "Map does not have an associated linear map", *this);
  613.   }
  614.  
  615. // Start to build matrix:
  616.  
  617.   Matrix hold = domain.EmbeddingSpace().HoistTangent();
  618.   Matrix holdr = range.EmbeddingSpace().HoistTangent();
  619.   Matrix alt = hold * t;
  620.   Matrix alinvt;      
  621.   if (invertible) alinvt = invt * Transpose(hold);
  622.  
  623. // The range depends on the character of the affine map.  If it is
  624. // to an affine subset of an affine space, the range is the
  625. // corresponding subspace of the tangent space.  If it is to an
  626. // affine subset of a vector space, the range is a vector
  627. // subspace in the vector space.  If it is to a vector subset,
  628. // this range is the same vector subset.
  629.  
  630.   SubSet alr;
  631.   GeObType alrt;
  632.  
  633.   if (range.EmbeddingSpace().Holds() == AFF_SPACE) {
  634.     alr = range.TangentSub();
  635.     alrt = AFF_VECTOR;
  636.     alt = alt * Transpose(holdr);
  637.     if (invertible) alinvt = holdr * alinvt;
  638.   } else if (range.Holds() == AFFINE_SUBSET) {
  639.     alr = range.TangentSub();
  640.     alrt = VECTOR;
  641.   } else {
  642.     alr = range;
  643.     alrt = VECTOR;
  644.   }
  645.  
  646.   Map retval(LIN_MAP, invertible, alr, domain.TangentSub(), alt,
  647.          alinvt, FALSE, FALSE, AFF_VECTOR, alrt);
  648.   return retval;  
  649. }
  650.  
  651. // ***********************************************************************
  652. //
  653. // If this map is an affine functional (i.e. an affine map into the reals),
  654. // this returns the vector in the dual of the tangent space corresponding
  655. // to the associated linear map.
  656. //
  657.  
  658. Vector Map::AssocDualVector(void)
  659. {
  660.   static char* sig = "Vector Map::AssocDualVector(void)";
  661.   Space res = range.EmbeddingSpace();
  662.  
  663.   if ((holds != AFF_MAP) || (res != Reals)) {
  664.     errh.ErrorExit(sig, "Not an affine map into the Reals", *this);
  665.   }
  666.  
  667.   Vector retval = Vector(this->AssocLinear());
  668.   return (retval);
  669. }
  670.  
  671.  
  672.  
  673. // ***********************************************************************
  674. // ***********************************************************************
  675. //
  676. // LinearMap Class
  677. //
  678. // ***********************************************************************
  679. // ***********************************************************************
  680. //
  681. // Public member functions
  682. //
  683. // ***********************************************************************
  684. //
  685. // A map can only be cast down to a Linear Map if it is one (no
  686. // conversions).
  687. //
  688.  
  689. LinearMap::LinearMap(Map &m) : (m)
  690. {
  691.   type = LIN_MAP;
  692.   if ((holds != LIN_MAP) && (holds != NULL_MAP)) {
  693.       errh.ErrorExit("LinearMap::LinearMap(Map &)",
  694.                      "Attempted to cast a non-linear map to a linear map", m);
  695.   }
  696. }
  697.  
  698. // ***********************************************************************
  699. //
  700. // Need to do memberwise assignment, since Matrix class members need to
  701. // do memberwise assignment:
  702. //
  703.  
  704. LinearMap& LinearMap::operator=(LinearMap& m)
  705. {
  706.   range = m.range;
  707.   domain = m.domain;
  708.   t = m.t;
  709.   invt = m.invt;
  710.   invertible = m.invertible;
  711.   holds = m.holds;
  712.   hasAL = m.hasAL;
  713.   isDefault = m.isDefault;
  714.   domtype = m.domtype;
  715.   rettype = m.rettype;
  716.   return (*this);
  717. }
  718.  
  719. // ***********************************************************************
  720. //
  721. // Output for debugging:
  722. //
  723.  
  724. void LinearMap::debug_out(ostream &c, int indent) 
  725. {
  726.   static char* name = "LinearMap Object\n";
  727.   this->data_out(c, indent, name);
  728.   return;
  729. }
  730.  
  731. // ***********************************************************************
  732. //
  733. // Simplest kind of map creation.  Since we are going between bases,
  734. // the map is invertible, and no subsets are involved:
  735. //
  736.  
  737. LinearMap::LinearMap(VBasis &b1, VBasis &b2)
  738. {
  739.   static char* sig = "LinearMap::LinearMap(VBasis&, VBasis&)";
  740.  
  741.   type = LIN_MAP;
  742.   holds = LIN_MAP;
  743.   domain = b1.SpaceOf().FullSet();
  744.   range = b2.SpaceOf().FullSet();
  745.   hasAL = FALSE;
  746.   isDefault = FALSE;
  747.   domtype = b1.SpaceOf().NativeType();
  748.   rettype = b2.SpaceOf().NativeType();
  749.  
  750.   if (range.Dim() != domain.Dim()) {
  751.     errh.ErrorExit(sig, "Range and domain do not have same dimension", b1, b2);
  752.   }
  753.  
  754.   t = b1.Fromstdb() * b2.Tostdb();
  755.   if (fabs(Det(t)) > EPSILON) {
  756.     invt = Inverse(t);
  757.     invertible = TRUE;
  758.   } else {
  759.     errh.ErrorExit(sig, "Unexpected error", b1, b2);
  760.   }
  761. }
  762.  
  763. // ***********************************************************************
  764. //
  765. // More general map creation.  The domain of the map is a whole space,
  766. // but the range can be a subset.  Each item in the GeObList represents
  767. // the image of the corresponding object in the basis.
  768.  
  769. LinearMap::LinearMap(VBasis& b, VSubSet& s, GeObList& v)
  770. {
  771.   static char* sig = "LinearMap::LinearMap(VBasis&, VSubSet&, GeObList&)";
  772.   Space rspace = s.EmbeddingSpace();
  773.  
  774.   type = LIN_MAP;
  775.   holds = LIN_MAP;
  776.   domain = b.SpaceOf().FullSet();
  777.   range = s;
  778.   hasAL = FALSE;
  779.   isDefault = FALSE;
  780.   domtype = b.SpaceOf().NativeType();
  781.   rettype = rspace.NativeType();
  782.  
  783. // Make sure that enough image objects have been specified:
  784.  
  785.   if (v.Length() != domain.Dim()) {
  786.     errh.ErrorExit(sig, "Incorrect number of objects specified", b, v);
  787.   }
  788.  
  789. // Map the objects into the range subset and use them to start building
  790. // the transform matrix:
  791.  
  792.   GeObType target = rspace.NativeType();
  793.   Matrix temp(v.Length(), rspace.Dim());
  794.  
  795.   if (!(this->ConvertList(v, target, range, &temp))) {
  796.     errh.ErrorExit(sig, "Object cannot be mapped into range subset", v, range);
  797.   } 
  798.  
  799. // Build the transform matrix, and then start to build the inverse by
  800. // converting the range objects to the subspace basis representation:
  801.  
  802.   t = b.Fromstdb() * temp;
  803.   temp = temp * range.FromFull();
  804.  
  805. // Now figure out if the map is invertible.  Two requirements must
  806. // be met.  The dimension of the range and domain subsets must match,
  807. // and the specified image objects in the range must be independent.
  808.  
  809.   if (range.Dim() != domain.Dim()) {
  810.     invertible = FALSE;
  811.   } else {
  812.     invertible = (fabs(Det(temp)) > EPSILON);
  813.   }
  814.  
  815.   if (invertible) {
  816.     invt = range.FromFull() * Inverse(temp) * b.Tostdb();
  817.   }
  818. }
  819.  
  820. // ***********************************************************************
  821. //
  822. // The reverse of the above constructor.
  823. //
  824.  
  825. LinearMap::LinearMap(VSubSet &s, GeObList &v, VBasis &b)
  826. {
  827.   static char* sig = "LinearMap::LinearMap(VSubSet&, GeObList&, VBasis&)";
  828.   Space dspace = s.EmbeddingSpace();
  829.  
  830.   type = LIN_MAP;
  831.   holds = LIN_MAP;
  832.   domain = s;
  833.   range = b.SpaceOf().FullSet();
  834.   hasAL = FALSE;
  835.   isDefault = FALSE;
  836.   domtype = dspace.NativeType();
  837.   rettype = b.SpaceOf().NativeType();
  838.  
  839. // Make sure that enough preimage objects have been specified:
  840.  
  841.   if (v.Length() != range.Dim()) {
  842.     errh.ErrorExit(sig, "Incorrect number of objects specified", b, v, s);
  843.   }
  844.   if (range.Dim() != domain.Dim()) {
  845.     errh.ErrorExit(sig, "Domain and range dimensions do not match", v, s, b);
  846.   }
  847.  
  848. // Map the objects into the domain subset and use them to start building
  849. // the transform matrix:
  850.  
  851.   GeObType target = dspace.NativeType();
  852.   Matrix temp1(v.Length(), dspace.Dim());
  853.  
  854.   if (!(this->ConvertList(v, target, domain, &temp1))) {
  855.     errh.ErrorExit(sig,
  856.            "Object cannot be mapped into domain subset", v, domain);
  857.   } 
  858.  
  859. // We now need to confirm that the preimage objects we have been given
  860. // span the domain subset.  If they do, build the transform matrix:
  861.  
  862.   Matrix temp2 = temp1 * domain.FromFull();
  863.   if (fabs(Det(temp2)) <= EPSILON) {
  864.     errh.ErrorExit(sig, "Objects in domain are not independent", v, s);
  865.   }
  866.   t = domain.FromFull() * Inverse(temp2) * b.Tostdb();
  867.  
  868. // Maps built this way must be invertible, so build the inverse:
  869.  
  870.   invertible = TRUE;
  871.   invt = b.Fromstdb() * temp1;
  872. }
  873.  
  874. // ***********************************************************************
  875. //
  876. // Here is the most general linear map constructor.
  877. //
  878.  
  879. LinearMap::LinearMap(VSubSet &s1, GeObList &v1, VSubSet &s2, GeObList &v2)
  880. {
  881.   static char* sig =
  882.     "LinearMap::LinearMap(VSubSet&, GeObList&, VSubSet&, GeObList&)";
  883.   Space dspace = s1.EmbeddingSpace();
  884.   Space rspace = s2.EmbeddingSpace();
  885.  
  886.   type = LIN_MAP;
  887.   holds = LIN_MAP;
  888.   domain = s1;
  889.   range = s2;
  890.   hasAL = FALSE;
  891.   isDefault = FALSE;
  892.   domtype = dspace.NativeType();
  893.   rettype = rspace.NativeType();
  894.  
  895. // Make sure that enough preimage and image objects have been specified:
  896.  
  897.   if (v1.Length() != domain.Dim()) {
  898.     errh.ErrorExit(sig, "Incorrect number of objects specified", s1, v1);
  899.   }
  900.   if (v2.Length() != v1.Length()) {
  901.     errh.ErrorExit(sig, "Incorrect number of objects specified", s1, s2, v2);
  902.   }
  903.  
  904. // Map the objects into the domain and range subsets and use them to start
  905. // building the transform matrix:
  906.  
  907.   GeObType target = dspace.NativeType();
  908.   Matrix tempd(v1.Length(), dspace.Dim());
  909.  
  910.   if (!(this->ConvertList(v1, target, domain, &tempd))) {
  911.     errh.ErrorExit(sig,
  912.            "Object cannot be mapped into domain subset", v1, domain);
  913.   } 
  914.  
  915.   target = rspace.NativeType();
  916.   Matrix tempr(v2.Length(), rspace.Dim());
  917.  
  918.   if (!(this->ConvertList(v2, target, range, &tempr))) {
  919.     errh.ErrorExit(sig,
  920.            "Object cannot be mapped into range subset", v2, range);
  921.   } 
  922.  
  923. // We now need to confirm that the preimage objects we have been given
  924. // span the domain subset.  If they do, build the transform matrix, and
  925. // then start to build the inverse by converting the range objects to
  926. // their subspace basis representation:
  927.  
  928.   Matrix tempd2 = tempd * domain.FromFull();
  929.   if (fabs(Det(tempd2)) <= EPSILON) {
  930.     errh.ErrorExit(sig, "Objects in domain are not independent", v1, s1);
  931.   }
  932.   t = domain.FromFull() * Inverse(tempd2) * tempr;
  933.   tempr = tempr * range.FromFull();
  934.  
  935. // Now figure out if the map is invertible.  Two requirements must
  936. // be met.  The dimension of the range and domain subsets must match,
  937. // and the specified image objects in the range must be independent.
  938.  
  939.   if (range.Dim() != domain.Dim()) {
  940.     invertible = FALSE;
  941.   } else {
  942.     invertible = (fabs(Det(tempr)) > EPSILON);
  943.   }
  944.  
  945.   if (invertible) {
  946.     invt = range.FromFull() * Inverse(tempr) * tempd;
  947.   }
  948. }
  949.  
  950. // ***********************************************************************
  951. // ***********************************************************************
  952. //
  953. // AffineMap Class
  954. //
  955. // ***********************************************************************
  956. // ***********************************************************************
  957. //
  958. // Private/protected member functions
  959. //
  960. // ***********************************************************************
  961. //
  962. // Builds the default standard affine maps between spaces in a space set.
  963. //
  964.  
  965. AffineMap::AffineMap(Space &s1, Space &s2)
  966. {
  967.   static char* sig = "AffineMap::AffineMap(Space&, Space&)";
  968.  
  969. // Standard affine maps go between standard affine subsets.  Determine
  970. // if these spaces have defined subsets; create them if they do not,
  971. // get them if they do.
  972.  
  973.   domain = s1.StdAffineSubset();
  974.   range = s2.StdAffineSubset(); 
  975.  
  976.   if (domain.Holds() == NULL_SUBSET) {
  977.     domain = ASubSet(s1);
  978.   }
  979.  
  980.   if (range.Holds() == NULL_SUBSET) {
  981.     range = ASubSet(s2);
  982.   }
  983.  
  984.   type = AFF_MAP;
  985.   holds = AFF_MAP;
  986.   invertible = TRUE;
  987.   isDefault = TRUE;
  988.   domtype = s1.NativeType();
  989.   rettype = s2.NativeType();
  990.  
  991.   SRel s1t = s1.ThisSpaceIs();
  992.   SRel s2t = s2.ThisSpaceIs();
  993.   Space check = s1.GetSpace(s2t);
  994.  
  995.   if (check != s2) {
  996.     errh.ErrorExit(sig, "Spaces are not linked", s1, s2);
  997.   }
  998.  
  999.   hasAL = (s1t == AFFINE) && (s2t == LINEARIZATION);
  1000.   t = IdentityMatrix(s1.MatrixSize());
  1001.   invt = t;
  1002. }
  1003.  
  1004. // ***********************************************************************
  1005. //
  1006. // Builds a consistent affine map for linking an affine space to
  1007. // a linearization space.  USER IS RESPONSIBLE FOR COPYING STD.
  1008. // AFFINE SUBSET GENERATED IN THIS ROUTINE TO THE EMBEDDING SPACE.
  1009.  
  1010.  
  1011. AffineMap::AffineMap(AffineMap& apm, ProjectiveMap& lpm)
  1012. {
  1013.   Space P = apm.Range().EmbeddingSpace();
  1014.   Space A = apm.Domain().EmbeddingSpace();
  1015.  
  1016. // Linearization space needs a consistent standard affine subset.
  1017.  
  1018.   range = A.StdAffineSubset();
  1019.   domain = ASubSet(P.StdAffineSubset(), lpm.Inv());
  1020.   type = AFF_MAP;
  1021.   holds = AFF_MAP;
  1022.   invertible = TRUE;
  1023.   hasAL = FALSE;
  1024.   isDefault = FALSE;
  1025.   domtype = VECTOR;
  1026.   rettype = AFF_POINT;
  1027.  
  1028.   t = lpm.t * apm.invt;
  1029.   invt = Inverse(t);
  1030. }
  1031.  
  1032. // ***********************************************************************
  1033. //
  1034. // Builds a consistent affine map for linking an affine space to
  1035. // a projective space:  
  1036. //
  1037.  
  1038. AffineMap::AffineMap(ProjectiveMap& lpm, AffineMap& lam)
  1039. {
  1040.   Space L = lpm.Domain().EmbeddingSpace();
  1041.   Space A = lam.Range().EmbeddingSpace();
  1042.  
  1043. // Projective space does not yet have a standard affine subset.  It
  1044. // must be consistent with the linearization space standard affine
  1045. // subset.  USER IS RESPONSIBLE FOR COPYING STD.
  1046. // AFFINE SUBSET GENERATED IN THIS ROUTINE TO THE EMBEDDING SPACE.
  1047.  
  1048.   domain = A.StdAffineSubset();
  1049.   range = ASubSet(L.StdAffineSubset(), lpm);
  1050.   type = AFF_MAP;
  1051.   holds = AFF_MAP;
  1052.   invertible = TRUE;
  1053.   hasAL = FALSE;
  1054.   isDefault = FALSE;
  1055.   domtype = AFF_POINT;
  1056.   rettype = PROJ_POINT;
  1057.  
  1058.   t = lam.invt * lpm.t;
  1059.   invt = Inverse(t);
  1060. }
  1061.  
  1062. // ***********************************************************************
  1063. //
  1064. // Public member functions
  1065. //
  1066. // ***********************************************************************
  1067. //
  1068. // A map can only be cast down to an affine map if it is one (no
  1069. // conversions). 
  1070. //
  1071.  
  1072. AffineMap::AffineMap(Map &m) : (m)
  1073. {
  1074.   type = AFF_MAP;
  1075.   if ((holds != AFF_MAP) && (holds != NULL_MAP)) {
  1076.       errh.ErrorExit("AffineMap::AffineMap(Map&)",
  1077.          "Attempted to cast a non-affine map to an affine map", m);
  1078.   }
  1079. }
  1080.  
  1081. // ***********************************************************************
  1082. //
  1083. // Need to do memberwise assignment, since Matrix class members need to
  1084. // do memberwise assignment.
  1085. //
  1086.  
  1087. AffineMap& AffineMap::operator=(AffineMap& m)
  1088. {
  1089.   range = m.range;
  1090.   domain = m.domain;
  1091.   t = m.t;
  1092.   invt = m.invt;
  1093.   invertible = m.invertible;
  1094.   holds = m.holds;
  1095.   hasAL = m.hasAL;
  1096.   isDefault = m.isDefault;
  1097.   domtype = m.domtype;
  1098.   rettype = m.rettype;
  1099.   return (*this);
  1100. }
  1101.  
  1102. // ***********************************************************************
  1103. //
  1104. // Output for debugging:
  1105. //
  1106.  
  1107. void AffineMap::debug_out(ostream &c, int indent) 
  1108. {
  1109.   static char* name = "AffineMap Object\n";
  1110.   this->data_out(c, indent, name);
  1111.   return;
  1112. }
  1113.  
  1114. // ***********************************************************************
  1115. //
  1116. // Simplest kind of map creation.  Since we are going between bases,
  1117. // the map is invertible, and no subsets are involved.  Restrict this
  1118. // constructor to just handle maps between affine spaces (i.e. no
  1119. // maps from simplex to vbasis).
  1120. //
  1121.  
  1122. AffineMap::AffineMap(Basis &b1, Basis &b2)
  1123. {
  1124.   static char* sig = "AffineMap::AffineMap(Basis&, Basis&)";
  1125.  
  1126.   type = AFF_MAP;
  1127.   holds = AFF_MAP;
  1128.   domain = b1.SpaceOf().FullSet();
  1129.   range = b2.SpaceOf().FullSet();
  1130.   hasAL = TRUE;
  1131.   isDefault = FALSE;
  1132.   domtype = AFF_POINT;
  1133.   rettype = AFF_POINT;
  1134.  
  1135. // Check that the type of basis matches:
  1136.  
  1137.   BasisType b1t = b1.Holds();
  1138.   BasisType b2t = b2.Holds();
  1139.  
  1140.   if ((b1t != b2t) || !((b1t == SIMPLEX) || (b1t == FRAME))) {
  1141.     errh.ErrorExit(sig, "Specified bases are not of the correct type", b1, b2);
  1142.   }
  1143.  
  1144. // Check that the dimensionality matches:
  1145.  
  1146.   if (range.Dim() != domain.Dim()) {
  1147.     errh.ErrorExit(sig, "Range and domain do not have same dimension", b1, b2);
  1148.   }
  1149.  
  1150.   t = b1.Fromstdb() * b2.Tostdb();
  1151.   if (fabs(Det(t)) > EPSILON) {
  1152.     invt = Inverse(t);
  1153.     invertible = TRUE;
  1154.   } else {
  1155.     errh.ErrorExit(sig, "Unexpected error", b1, b2);
  1156.   }
  1157. }
  1158.  
  1159. // ***********************************************************************
  1160. //
  1161. // Map from a full affine space into some affine or vector subset.
  1162. // We currently require that a simplex be used to define the domain
  1163. // preimage.
  1164. //
  1165.  
  1166. AffineMap::AffineMap(Simplex& b, SubSet& s, GeObList& v)
  1167. {
  1168.   static char* sig = "AffineMap::AffineMap(Simplex&, SubSet&, GeObList&)";
  1169.   Space rspace = s.EmbeddingSpace();
  1170.  
  1171.   type = AFF_MAP;
  1172.   holds = AFF_MAP;
  1173.   domain = b.SpaceOf().FullSet();
  1174.   range = s;
  1175.   domtype = AFF_POINT;
  1176.   rettype = rspace.NativeType();
  1177.   hasAL = (rettype != PROJ_POINT);
  1178.   isDefault = FALSE;
  1179.  
  1180. // Make sure the subset is an affine or linear subset and that enough
  1181. // image objects have been specified:
  1182.  
  1183.   if ((s.Holds() != LINEAR_SUBSET) && (s.Holds() != AFFINE_SUBSET)) {
  1184.     errh.ErrorExit(sig, "Specified subset is not linear or affine", s);
  1185.   }
  1186.   if (v.Length() != domain.Dim() + 1) {
  1187.     errh.ErrorExit(sig, "Incorrect number of objects specified", b, v);
  1188.   }
  1189.  
  1190. // Map the objects into the range subset and use them to start building
  1191. // the transform matrix:
  1192.  
  1193.   GeObType target = rspace.NativeType();
  1194.   Matrix temp(v.Length(), rspace.MatrixSize());
  1195.  
  1196.   if (!(this->ConvertList(v, target, range, &temp))) {
  1197.     errh.ErrorExit(sig, "Object cannot be mapped into range subset", v, range);
  1198.   } 
  1199.  
  1200. // Build the transform matrix, and then start to build the inverse by
  1201. // converting the range objects to the subspace basis representation:
  1202.  
  1203.   t = b.Fromstdb() * temp;
  1204.   temp = temp * range.FromFull();
  1205.  
  1206. // Now figure out if the map is invertible.  Three requirements must
  1207. // be met.  The dimension of the range and domain subsets must match,
  1208. // the range must be an affine subset, and the specified image objects
  1209. // in the range must be independent.
  1210.  
  1211.   if ((range.Dim() != domain.Dim()) || (range.Holds() != AFFINE_SUBSET)) {
  1212.     invertible = FALSE;
  1213.   } else {
  1214.     invertible = (fabs(Det(temp)) > EPSILON);
  1215.   }
  1216.  
  1217.   if (invertible) {
  1218.     invt = range.FromFull() * Inverse(temp) * b.Tostdb();
  1219.   }
  1220. }
  1221.  
  1222. // ***********************************************************************
  1223. //
  1224. // Map from an affine subset to a full affine space. 
  1225. //
  1226.  
  1227. AffineMap::AffineMap(ASubSet& s, GeObList& v, Simplex& b)
  1228. {
  1229.   static char* sig = "AffineMap::AffineMap(ASubSet&, GeObList&, Simplex&)";
  1230.   Space dspace = s.EmbeddingSpace();
  1231.  
  1232.   type = AFF_MAP;
  1233.   holds = AFF_MAP;
  1234.   domain = s;
  1235.   range = b.SpaceOf().FullSet();
  1236.   domtype = dspace.NativeType();
  1237.   rettype = AFF_POINT;
  1238.   hasAL = (domtype == AFF_POINT);
  1239.   isDefault = FALSE;
  1240.  
  1241. // Make sure the subset is an affine subset and that enough
  1242. // preimage objects have been specified:
  1243.  
  1244.   if (v.Length() != range.Dim() + 1) {
  1245.     errh.ErrorExit(sig, "Incorrect number of objects specified", b, v, s);
  1246.   }
  1247.   if (range.Dim() != domain.Dim()) {
  1248.     errh.ErrorExit(sig, "Domain and range dimensions do not match", v, s, b);
  1249.   }
  1250.  
  1251. // Map the objects into the domain subset and use them to start building
  1252. // the transform matrix:
  1253.  
  1254.   GeObType target = dspace.NativeType();
  1255.   Matrix temp1(v.Length(), dspace.MatrixSize());
  1256.  
  1257.   if (!(this->ConvertList(v, target, domain, &temp1))) {
  1258.     errh.ErrorExit(sig,
  1259.            "Object cannot be mapped into domain subset", v, domain);
  1260.   } 
  1261.  
  1262. // We now need to confirm that the preimage objects we have been given
  1263. // span the domain subset.  If they do, build the transform matrix:
  1264.  
  1265.   Matrix temp2 = temp1 * domain.FromFull();
  1266.   if (fabs(Det(temp2)) <= EPSILON) {
  1267.     errh.ErrorExit(sig, "Objects in domain are not independent", v, s);
  1268.   }
  1269.   t = domain.FromFull() * Inverse(temp2) * b.Tostdb();
  1270.  
  1271. // Maps built this way must be invertible, so build the inverse:
  1272.  
  1273.   invertible = TRUE;
  1274.   invt = b.Fromstdb() * temp1;
  1275. }
  1276.  
  1277. // ***********************************************************************
  1278. //
  1279. // Most general constructor for affine maps. 
  1280. //
  1281.  
  1282. AffineMap::AffineMap(ASubSet &s1, GeObList &v1, SubSet &s2, GeObList &v2)
  1283. {
  1284.   static char* sig =
  1285.     "AffineMap::AffineMap(ASubSet&, GeObList&, SubSet&, GeObList&)";
  1286.   Space dspace = s1.EmbeddingSpace();
  1287.   Space rspace = s2.EmbeddingSpace();
  1288.  
  1289.   type = AFF_MAP;
  1290.   holds = AFF_MAP;
  1291.   domain = s1;
  1292.   range = s2;
  1293.  
  1294.   domtype = dspace.NativeType();
  1295.   rettype = rspace.NativeType();
  1296.   hasAL = ((domtype == AFF_POINT) && (rettype != PROJ_POINT));
  1297.   isDefault = FALSE;
  1298.  
  1299. // The range subspace can be either affine or linear.  Also make
  1300. // sure enough objects have been specified in each list:
  1301.  
  1302.   if ((s2.Holds() != LINEAR_SUBSET) && (s2.Holds() != AFFINE_SUBSET)) {
  1303.     errh.ErrorExit(sig, "Range subset is not linear or affine", s2);
  1304.   }
  1305.   if (v1.Length() != domain.Dim() + 1) {
  1306.     errh.ErrorExit(sig, "Incorrect number of objects specified", s1, v1);
  1307.   }
  1308.   if (v2.Length() != v1.Length()) {
  1309.     errh.ErrorExit(sig, "Incorrect number of objects specified", v1, s2, v2);
  1310.   }
  1311.  
  1312. // Map the objects into the domain and range subsets and use them to start
  1313. // building the transform matrix:
  1314.  
  1315.   GeObType target = dspace.NativeType();
  1316.   Matrix tempd(v1.Length(), dspace.MatrixSize());
  1317.  
  1318.   if (!(this->ConvertList(v1, target, domain, &tempd))) {
  1319.     errh.ErrorExit(sig,
  1320.            "Object cannot be mapped into domain subset", v1, domain);
  1321.   } 
  1322.  
  1323.   target = rspace.NativeType();
  1324.   Matrix tempr(v2.Length(), rspace.MatrixSize());
  1325.  
  1326.   if (!(this->ConvertList(v2, target, range, &tempr))) {
  1327.     errh.ErrorExit(sig,
  1328.            "Object cannot be mapped into range subset", v2, range);
  1329.   } 
  1330.  
  1331. // We now need to confirm that the preimage objects we have been given
  1332. // span the domain subset.  If they do, build the transform matrix, and
  1333. // then start to build the inverse by converting the range objects to
  1334. // their subspace basis representation:
  1335.  
  1336.   Matrix tempd2 = tempd * domain.FromFull();
  1337.   if (fabs(Det(tempd2)) <= EPSILON) {
  1338.     errh.ErrorExit(sig, "Objects in domain are not independent", v1, s1);
  1339.   }
  1340.   t = domain.FromFull() * Inverse(tempd2) * tempr;
  1341.   tempr = tempr * range.FromFull();
  1342.  
  1343. // Now figure out if the map is invertible.  Two requirements must
  1344. // be met.  The dimension of the range and domain subsets must match,
  1345. // and the specified image objects in the range must be independent.
  1346.  
  1347.   if (range.Dim() != domain.Dim()) {
  1348.     invertible = FALSE;
  1349.   } else {
  1350.     invertible = (fabs(Det(tempr)) > EPSILON);
  1351.   }
  1352.  
  1353.   if (invertible) {
  1354.     invt = range.FromFull() * Inverse(tempr) * tempd;
  1355.   }
  1356. }
  1357.  
  1358. // ***********************************************************************
  1359. // ***********************************************************************
  1360. //
  1361. // ProjectiveMap Class
  1362. //
  1363. // ***********************************************************************
  1364. // ***********************************************************************
  1365. //
  1366. // Private/protected member functions
  1367. //
  1368. // ***********************************************************************
  1369. //
  1370. // A core problem when constructing projective maps is scaling the
  1371. // rows of the transformation matrix so the "extra" point maps
  1372. // to its image.  This routine scales the rows of a matrix so the
  1373. // specified "unit point" matrix is the sum of the rows.  It does this
  1374. // by solving a system of equations:
  1375. //
  1376.  
  1377. Boolean ProjectiveMap::ScaleRows(Matrix* mat, Matrix& unit_pt)
  1378. {
  1379.   if (fabs(Det(*mat)) < EPSILON) {
  1380.     return (FALSE);
  1381.   }
  1382.   Matrix invmat = Inverse(*mat);
  1383.   
  1384. // Solve for the row coefficients.  If any coefficient = 0, the
  1385. // unit point was not in general position.  Otherwise, multiply
  1386. // the matrix rows by the coefficients, and report the success:
  1387.  
  1388.   Matrix coeff = unit_pt * invmat;
  1389.   for (int i = 0; i < coeff.Columns(); i++) {
  1390.     if (fabs(coeff[0][i]) < EPSILON) {
  1391.       return (FALSE);
  1392.     } else {
  1393.       (*mat)[i] = (coeff[0][i] * (*mat)[i])[0];
  1394.     }
  1395.   }
  1396.   return (TRUE);
  1397. }
  1398.  
  1399. // ***********************************************************************
  1400. //
  1401. // Similar to the general ConvertList routine, but sticks the last
  1402. // object into a separate Matrix:
  1403. //
  1404.  
  1405. Boolean ProjectiveMap::ConvertListUnit(GeObList& v, GeObType trg, SubSet& d,
  1406.                        Matrix* temp, Matrix* unit)
  1407. {
  1408.   Space destspace = d.EmbeddingSpace();
  1409.  
  1410.   for (int i = 0; i < v.Length(); i++) {
  1411.     GeOb hold = v[i].MapTo(trg);
  1412.     if (hold.SpaceOf() != destspace) {
  1413.       return (FALSE);
  1414.     } 
  1415.     if (!d.IsIn(hold)) {
  1416.       return (FALSE);
  1417.     }
  1418.     if (i == v.Length() - 1) {
  1419.       (*unit)[0] = hold.MatrixRep()[0];
  1420.     } else {
  1421.       (*temp)[i] = hold.MatrixRep()[0];
  1422.     }
  1423.   }
  1424.   return (TRUE);
  1425. }
  1426.  
  1427.  
  1428. // ***********************************************************************
  1429. //
  1430. // Builds the default standard projective maps between spaces in a space set.
  1431. // One space must be projective, the other a linearization space.
  1432.  
  1433. ProjectiveMap::ProjectiveMap(Space& s1, Space& s2)
  1434. {
  1435.   static char* sig = "ProjectiveMap::ProjectiveMap(Space&, Space&)";
  1436.  
  1437.   domain = s1.FullProjSet();
  1438.   range = s2.FullProjSet();
  1439.  
  1440.   type = PROJ_MAP;
  1441.   holds = PROJ_MAP;
  1442.   invertible = TRUE;
  1443.   hasAL = FALSE;
  1444.   isDefault = TRUE;
  1445.   domtype = domain.Accepts();
  1446.   rettype = range.Accepts();
  1447.  
  1448.   SRel s1t = s1.ThisSpaceIs();
  1449.   SRel s2t = s2.ThisSpaceIs();
  1450.   Space check = s1.GetSpace(s2t);
  1451.  
  1452.   if (check != s2) {
  1453.     errh.ErrorExit(sig, "Spaces are not linked", s1, s2);
  1454.   }
  1455.  
  1456.   t = IdentityMatrix(s1.MatrixSize());
  1457.   invt = t;
  1458. }
  1459.  
  1460. // ***********************************************************************
  1461. //
  1462. // Builds a consistent projective map for linking a projective space to
  1463. // a linearization space.
  1464. //
  1465.  
  1466. ProjectiveMap::ProjectiveMap(VSpace& L, AffineMap& lpm, PSpace& P)
  1467. {
  1468.   domain = L.FullProjSet();
  1469.   range = P.FullSet();
  1470.   type = PROJ_MAP;
  1471.   holds = PROJ_MAP;
  1472.   invertible = TRUE;
  1473.   hasAL = FALSE;
  1474.   isDefault = lpm.isDefault;
  1475.   domtype = VEC_EC;
  1476.   rettype = PROJ_POINT;
  1477.  
  1478.   t = lpm.t;
  1479.   invt = Inverse(t);
  1480. }
  1481.  
  1482. // ***********************************************************************
  1483. //
  1484. // Public member functions
  1485. //
  1486. // ***********************************************************************
  1487. //
  1488. // A map can only be cast down to a Projective Map if it is one (no
  1489. // conversions). 
  1490.  
  1491. ProjectiveMap::ProjectiveMap(Map &m) : (m)
  1492. {
  1493.   type = PROJ_MAP;
  1494.   if ((holds != PROJ_MAP) && (holds != NULL_MAP)) {
  1495.       errh.ErrorExit("ProjectiveMap::ProjectiveMap(Map&)",
  1496.          "Attempted to cast a non-projective map to a projective map", m);
  1497.   }
  1498. }
  1499.  
  1500. // ***********************************************************************
  1501. //
  1502. // Need to do memberwise assignment, since Matrix class members need to
  1503. // do memberwise assignment:
  1504. //
  1505.  
  1506. ProjectiveMap& ProjectiveMap::operator=(ProjectiveMap& m)
  1507. {
  1508.   range = m.range;
  1509.   domain = m.domain;
  1510.   t = m.t;
  1511.   invt = m.invt;
  1512.   invertible = m.invertible;
  1513.   holds = m.holds;
  1514.   hasAL = m.hasAL;
  1515.   isDefault = m.isDefault;
  1516.   domtype = m.domtype;
  1517.   rettype = m.rettype;
  1518.   return (*this);
  1519. }
  1520.  
  1521. // ***********************************************************************
  1522. //
  1523. // Debug output
  1524. //
  1525.  
  1526. void ProjectiveMap::debug_out(ostream &c, int indent) 
  1527. {
  1528.   static char* name = "ProjectiveMap Object\n";
  1529.   this->data_out(c, indent, name);
  1530.   return;
  1531. }
  1532.  
  1533. // ***********************************************************************
  1534. //
  1535. // Used for invertible maps between whole projective spaces:
  1536. //
  1537.  
  1538. ProjectiveMap::ProjectiveMap(HFrame &b1, HFrame &b2)
  1539. {
  1540.   static char* sig = "ProjectiveMap::ProjectiveMap(HFrame&, HFrame&)";
  1541.  
  1542.   type = PROJ_MAP;
  1543.   holds = PROJ_MAP;
  1544.   domain = b1.SpaceOf().FullSet();
  1545.   range = b2.SpaceOf().FullSet();
  1546.   hasAL = FALSE;
  1547.   isDefault = FALSE;
  1548.   domtype = PROJ_POINT;
  1549.   rettype = PROJ_POINT;
  1550.  
  1551. // Check that the type of basis matches:
  1552.  
  1553.   if ((b1.Holds() != HFRAME) || (b2.Holds() != HFRAME)) {
  1554.     errh.ErrorExit(sig, "Specified bases are not projective frames", b1, b2);
  1555.   }
  1556.  
  1557. // Check that the dimensionality matches:
  1558.  
  1559.   if (range.Dim() != domain.Dim()) {
  1560.     errh.ErrorExit(sig, "Range and domain do not have same dimension", b1, b2);
  1561.   }
  1562.  
  1563.   t = b1.Fromstdb() * b2.Tostdb();
  1564.   if (fabs(Det(t)) > EPSILON) {
  1565.     invt = Inverse(t);
  1566.     invertible = TRUE;
  1567.   } else {
  1568.     errh.ErrorExit(sig, "Unexpected error", b1, b2);
  1569.   }
  1570. }
  1571.  
  1572. // ***********************************************************************
  1573. //
  1574. // Used for invertible maps from a whole projective space to
  1575. // a projective subset.  Since the only allowable method for
  1576. // creating non-invertible maps is to have a domain subset with
  1577. // removed points, we require the image objects to be independent.
  1578. // Note that the subset could be from a vector space.  The range
  1579. // subset cannot have removed points, since there would be no way
  1580. // of returning a unique point in the embedding space.
  1581. //
  1582.  
  1583. ProjectiveMap::ProjectiveMap(HFrame &b, PSubSet &s, GeObList &v)
  1584. {
  1585.   static char* sig =
  1586.       "ProjectiveMap::ProjectiveMap(HFrame&, PSubSet&, GeObList&)";
  1587.   Space rspace = s.EmbeddingSpace();
  1588.  
  1589.   type = PROJ_MAP;
  1590.   holds = PROJ_MAP;
  1591.   domain = b.SpaceOf().FullSet();
  1592.   range = s;
  1593.   hasAL = FALSE;
  1594.   isDefault = FALSE;
  1595.   domtype = PROJ_POINT;
  1596.   rettype = s.Accepts();
  1597.  
  1598. // Make that enough image objects have been specified.  Also, since 
  1599. // only invertible maps can be specified, check the dimensions:
  1600.  
  1601.   if (v.Length() != domain.Dim() + 2) {
  1602.     errh.ErrorExit(sig, "Incorrect number of objects specified", b, v);
  1603.   }
  1604.   if (range.Dim() != domain.Dim()) {
  1605.     errh.ErrorExit(sig, "Range and domain dimensions must be equal",
  1606.             range, domain);
  1607.   }
  1608.   if (range.HasRemovedPoints()) {
  1609.     errh.ErrorExit(sig, "Range subset cannot have removed points", range);
  1610.   }
  1611.  
  1612. // Map the objects into the range subset and use them to start building
  1613. // the transform matrix:
  1614.  
  1615.   int rdim = rspace.MatrixSize(); 
  1616.   Matrix temp(v.Length() - 1, rdim);
  1617.   Matrix unit(1, rdim);
  1618.  
  1619.   if (!(this->ConvertListUnit(v, rettype, range, &temp, &unit))) {
  1620.     errh.ErrorExit(sig, "Object cannot be mapped into range subset", v, range);
  1621.   } 
  1622.  
  1623. // Convert the range objects to the range subset basis representation, then
  1624. // scale the matrix rows:
  1625.  
  1626.   temp = temp * range.FromFull();
  1627.   unit = unit * range.FromFull();
  1628.   if (!(this->ScaleRows(&temp, unit))) {
  1629.     errh.ErrorExit(sig, "Range objects not in general position", v, s);
  1630.   }
  1631.  
  1632. // Build the transform matrix, then build the inverse:
  1633.  
  1634.   t = b.Fromstdb() * temp * range.ToFull();
  1635.   invt = range.FromFull() * Inverse(temp) * b.Tostdb();
  1636.   invertible = TRUE;
  1637. }
  1638.  
  1639. // ***********************************************************************
  1640. //
  1641. // Used for map from a projective subset to a projective space.  If
  1642. // subset has removed points, map will be non-invertible:
  1643. //
  1644.  
  1645. ProjectiveMap::ProjectiveMap(PSubSet &s, GeObList &v, HFrame &b)
  1646. {
  1647.   static char* sig =
  1648.      "ProjectiveMap::ProjectiveMap(PSubSet&, GeObList&, HFrame&)";
  1649.   Space dspace = s.EmbeddingSpace();
  1650.  
  1651.   type = PROJ_MAP;
  1652.   holds = PROJ_MAP;
  1653.   domain = s;
  1654.   range = b.SpaceOf().FullSet();
  1655.   hasAL = FALSE;
  1656.   isDefault = FALSE;
  1657.   domtype = s.Accepts();
  1658.   rettype = PROJ_POINT;
  1659.  
  1660. // Make sure that enough preimage objects have been specified.
  1661. // Check that the domain and range dimensions match.
  1662.  
  1663.   if (v.Length() != range.Dim() + 2) {
  1664.     errh.ErrorExit(sig, "Incorrect number of objects specified", b, v, s);
  1665.   }
  1666.   if (range.Dim() != domain.Dim()) {
  1667.     errh.ErrorExit(sig, "Domain and range dimensions do not match", v, s, b);
  1668.   }
  1669.  
  1670. // Map the objects into the domain subset and use them to start building
  1671. // the transform matrix:
  1672.  
  1673.   int ddim = dspace.MatrixSize();
  1674.   Matrix temp(v.Length() - 1, ddim);
  1675.   Matrix unit(1, ddim);
  1676.  
  1677.   if (!(this->ConvertListUnit(v, domtype, domain, &temp, &unit))) {
  1678.     errh.ErrorExit(sig,
  1679.            "Object cannot be mapped into domain subset", v, domain);
  1680.   } 
  1681.  
  1682. // Convert the domain objects to the domain subset basis representation, then
  1683. // scale the matrix rows:
  1684.  
  1685.   temp = temp * domain.FromFull();
  1686.   unit = unit * domain.FromFull();
  1687.   if (!(this->ScaleRows(&temp, unit))) {
  1688.     errh.ErrorExit(sig, "Domain objects not in general position", v, s);
  1689.   } 
  1690.  
  1691. // Build the transform matrix:
  1692.  
  1693.   t = domain.FromFull() * Inverse(temp) * b.Tostdb();
  1694.  
  1695. // The invertibility of this transform depends on whether the domain
  1696. // subset has been defined to have removed points.  If it has, the
  1697. // map is not invertible.  Otherwise, calculate the inverse:
  1698.  
  1699.   invertible = !domain.HasRemovedPoints();
  1700.  
  1701.   if (invertible) {
  1702.     invt = b.Fromstdb() * temp * domain.ToFull();
  1703.   }
  1704. }
  1705.  
  1706. // ***********************************************************************
  1707. //
  1708. // Used for maps between subsets.  If domain subset has removed points,
  1709. // map will not be invertible:
  1710. //
  1711.  
  1712. ProjectiveMap::ProjectiveMap(PSubSet &s1, GeObList &v1,
  1713.                  PSubSet &s2, GeObList &v2)
  1714. {
  1715.   static char* sig =
  1716.       "ProjectiveMap::ProjectiveMap(PSubSet&, GeObList&, PSubSet&, GeObList&)";
  1717.   Space dspace = s1.EmbeddingSpace();
  1718.   Space rspace = s2.EmbeddingSpace();
  1719.  
  1720.   type = PROJ_MAP;
  1721.   holds = PROJ_MAP;
  1722.   domain = s1;
  1723.   range = s2;
  1724.   hasAL = FALSE;
  1725.   isDefault = FALSE;
  1726.   domtype = s1.Accepts();
  1727.   rettype = s2.Accepts();
  1728.  
  1729. // Make sure that enough preimage and image objects have been specified.
  1730. // Check that the domain and range dimensions match.
  1731.  
  1732.   if (v1.Length() != domain.Dim() + 2) {
  1733.     errh.ErrorExit(sig, "Incorrect number of objects specified", v1, s1);
  1734.   }
  1735.   if (v2.Length() != range.Dim() + 2) {
  1736.     errh.ErrorExit(sig, "Incorrect number of objects specified", v2, s2);
  1737.   }
  1738.   if (range.Dim() != domain.Dim()) {
  1739.     errh.ErrorExit(sig, "Domain and range dimensions do not match", s1, s2);
  1740.   }
  1741.   if (range.HasRemovedPoints()) {
  1742.     errh.ErrorExit(sig, "Range subset cannot have removed points", range);
  1743.   }
  1744.  
  1745. // Map the objects into the domain and range subsets and use them to start
  1746. // building the transform matrix:
  1747.  
  1748.   int ddim = dspace.MatrixSize();
  1749.   Matrix tempd(v1.Length() - 1, ddim);
  1750.   Matrix unitd(1, ddim);
  1751.   if (!(this->ConvertListUnit(v1, domtype, domain, &tempd, &unitd))) {
  1752.     errh.ErrorExit(sig,
  1753.            "Object cannot be mapped into domain subset", v1, domain);
  1754.   } 
  1755.  
  1756.   int rdim = rspace.MatrixSize();
  1757.   Matrix tempr(v2.Length() - 1, rdim);
  1758.   Matrix unitr(1, rdim);
  1759.   if (!(this->ConvertListUnit(v2, rettype, range, &tempr, &unitr))) {
  1760.     errh.ErrorExit(sig,
  1761.            "Object cannot be mapped into range subset", v2, range);
  1762.   }
  1763.  
  1764. // Convert the domain objects to the domain subset basis representation, then
  1765. // scale the matrix rows.  Do the same for the range:
  1766.  
  1767.   tempd = tempd * domain.FromFull();
  1768.   unitd = unitd * domain.FromFull();
  1769.   if (!(this->ScaleRows(&tempd, unitd))) {
  1770.     errh.ErrorExit(sig, "Domain objects not in general position", v1, s1);
  1771.   } 
  1772.  
  1773.   tempr = tempr * range.FromFull();
  1774.   unitr = unitr * range.FromFull();
  1775.   if (!(this->ScaleRows(&tempr, unitr))) {
  1776.     errh.ErrorExit(sig, "Range objects not in general position", v2, s2);
  1777.   } 
  1778.  
  1779. // Build the transform matrix:
  1780.  
  1781.   t = domain.FromFull() * Inverse(tempd) * tempr * range.ToFull();
  1782.  
  1783. // The invertibility of this transform depends on whether the domain
  1784. // subset has been defined to have removed points.  If it has, the
  1785. // map is not invertible.  Otherwise, calculate the inverse:
  1786.  
  1787.   invertible = !domain.HasRemovedPoints();
  1788.  
  1789.   if (invertible) {
  1790.     invt = range.FromFull() * Inverse(tempr) * tempd * domain.ToFull();
  1791.   }
  1792. }
  1793. // ***********************************************************************
  1794.